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ABSTRACT 

Context. Embedded planets disturb the density structure of the ambient disk, and gravitational back-reaction possibly will induce a 
change in the planet's orbital elements. Low-mass planets only have a weak impact on the disk, so their wake's torque can be treated 
in linear theory. Larger planets will begin to open up a gap in the disk through nonlinear interaction. Accurate determination of the 
l/") , forces acting on the planet requires careful numerical analysis. Recently, the validity of the often used fast orbital advection algorithm 

(FARGO) has been put into question, and special numerical resolution and stability requirements have been suggested. 
Aims. We study the process of planet-disk interaction for low-mass planets of a few Earth masses, and reanalyze the numerical 
requirements to obtain converged and stable results. One focus lies on the applicability of the FARGO-algorithm. Additionally, we 
study the difference of two and three-dimensional simulations, compare global with local setups, as well as isothermal and adiabatic 
conditions. 

Methods. We study the influence of the planet on the disk through two- and three-dimensional hydrodynamical simulations. To 
strengthen our conclusions we perform a detailed numerical comparison where several upwind and Riemann-solver based codes are 
used with and without the FARGO-algorithm. 

Results. With respect to the wake structure and the torque density acting on the planet, we demonstrate that the FARGO-algorithm 
?— j ' yields correct a correct and stable evolution for the planet-disk problem, and that at a fraction of the regular cpu-time. We find that 

C/j ' the resolution requirements for achieving convergent results in unshocked regions are rather modest and depend on the pressure scale 

height H of the disk. By comparing the torque densities of two- and three-dimensional simulations we show that a suitable vertical 
averaging procedure for the force gives an excellent agreement between the two. We show that isothermal and adiabatic runs can 
differ considerably, even for adiabatic indices very close to unity. 

^ ' Key words, accretion, accretion disks - protoplanetary disks - hydrodynamics - methods: numerical - planets and satellites: forma- 

ts ' tion 

^] ', 1. Introduction varying timestep size caused by the differentially rotating disk. 

OO ' Because the disk is highly supersonic with (azimuthal) Mach 

Q ■ Very young planets that are still embedded in the protoplanetary num bers of about 10 to 50, the angular velocity at the inner 

(S| ; disk will disturb the ambient density by their gravity. This will disk win i imit the timestep of the whole simulation, even though 

. lead to gravitational torques that can alter the orbital elements the planet or other reg ions of interest are located much farther 

L; ■ of the planet. For massive enough planets, the wake becomes out Changing to a rotating coordinate system will not help too 

. ,h ! nonlinear, and gap formation sets in. In numerical simulations much owing t0 the strong differential shear. T o solve t his par- 

^ . of embedded planets, different length scales have to be resolved, ticular pro blem and speed up the computation, iMassetl (l2000ih 

»H ; in particular when studying low-mass planets. On the one hand, has developed a fast orbital advection algorithm (FARGO). This 

, P. , the global structure has to be resolved to be able to obtain the met hod consists of an analytic, exact shift in the hydrodynami- 

correct structure of the wakes, i.e. the spiral arms generated by cal quan tities by approximately the average azimuthal velocity, 

the planet, which requires a sufficiently large radial domain. The The transport step utilizes only the residual velocity, which is 

libration of co-orbital material on horseshoe streamlines requires c i ose to the local sound speed. Depending on the grid layout and 

a full azimuthal extent of 2 n radians to be properly captured. On t he chosen radial range, a very large speed-up can be achieved, 

the other hand, the direct vicinity of the planet has to be resolved while at me same time th e intrinsic num erical diffusion of the 

to study detail effects, such as horseshoe drag or accretion onto sc heme is highly reduced ( lMassetll2000allbl) . 
the planet. 

To ease computational requirements, often planet-disk sim- The original version of the algorithm has been implemented 
ulations are performed in the two-dimensional (2D) thin disk into the public code FARGO, which is very often used in planet- 
approximation, because a full three-dimensional (3D) treatment disk and related simulations. The accuracy of the FARGO- 
with high resolution is still very time-consuming. However, even algorithm has been demonstrated in a detailed planet-disk com- 
under this reduced dimensionality, the problem is still com- parison project utilizing embedde d Neptune and Jupiter mass 
putationally very demanding. The main reason is the strongly planets dde Val-Borro et ail 120061) . There, it has been shown 
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that it leads to identical density profiles near the planet and 
total torques acting on the planet. Meanwhile, similar or- 
bital advection algorithms have been implemented into a va- 
riety of differ ent codes in two and three spatial dim ensions, 
e.g. NIRVANA dZiegler &Yorkell 19971: iKlev et al.ll2009l) . ATHENA 
([Gardiner & Stone! 1 120081: IStone & Gardiner! 1201 Ol). "and PLUTO 
dMignone et al.ll2007l |2012|) . Despite these widespread applica- 
tions, it has been claimed recently that usage of the FARGO- 
algorithm (here in connection with ATHENA) may lead to an 
unsteady behavior of the flow near the planet, which even af- 
fects the wake struc ture of the flow farther away from the planet 
dDong et al.H20TTbh . 

In the same paper, IDong et ail (1201 lbl) note that a very 
high numerical grid resolution is required to obtain a resolved 
flow near the planet. In particular, they analyze the smooth, un- 
shocked wake structure close to the planet and infer that a mini- 
mum spatial resolution of about 256 gridcells per scale height 
H, of the disk is needed to obtain good agreement with lin- 
ear studies. New simulations with a moving mesh technique 
also seem to indicate the n ecessity of very high resolutions 
dDuffell & MacFadvenll2012l) . 

Because a robust, fast, and reliable solution technique is 
mandatory in these type of simulations, we decided to address 
the planet-disk problem for a well defined standard setu p, which 
is very close to the one used in IDong et all (1201 lbl) . To an- 
swer the question of the validity of the FARGO-algorithm and 
estimate the resolution requirements, we applied several dif- 
ferent codes to an identical problem. These range from classi- 
cal second-order upwind schemes (e.g. RH2D, FARGO) to mod- 
ern Riemann-solvers such as PLUTO. The characteristics of these 
codes are specified in Appendix lAl 

Another critical issue in planet-disk simulations is the selec- 
tion of a realistic treatment of the gravitational force between the 
disk and the planet. Because the planet is typically treated as a 
pointmass and located within the numerical grid, regularization 
of the potential is required. In addition, physical smoothing is 
required to account for the otherwise neglected vertical thick- 
ness of the disk. The magnitude of this smoothing is hig hly rele- 
vant, since it influences the torques acting on the planet dMassetl 
120021) . and the smo othing parameter has even entered analyti- 
cal to rque formulas dMasset & Cas oli 2009:; iPaardekooper et all 
120101) . Because the 2D equations are obtained by a vertical av- 
eraging procedure, the force should be calculated by a suitable 
vertical inte gration as well. Thi s approach has been undertaken 
recently by iMiiller et alj d2012l) . who show that the smoothing 
length is indeed determined by the vertical thickness of the disk, 
and is roughly on the order of 0.7 H. They show in addition that 
the change of the disk thickness induced by the presence of the 
planet has to be taken into account. Because in recent simula- 
tions very short sm oothin g lengths have been used in 2 D simula- 
tions dDong et al.ll201 lbHDuffell & MacFadvenl2012l) , we com- 
pare our 2D results on the standard problem to an equivalent 3D 
setup and infer the required right amount of smoothing. 

Finally, we performed additional simulations for different 
equations of state. The first set of simulations deals with the of- 
ten used locally isothermal setup, while in comparison simula- 
tions we explore the outcome of adiabatic runs. This is impor- 
tant because some codes may not allow for treating an isother- 
mal equation of state. Here, we use different values for the ra- 
tio of specific heat y. In particular, a value of y very close to 
unity has often been quoted as closely resembling the isother- 
mal case. We show that this statement can depend on the physi- 
cal problem. In particular, in flows where the conservation of the 
entropy along streamlines is relevant, there can be strong differ- 



ences between an isothermal and an adiabatic flow, regardless 
of the value chosen fo r y. For the planet-disk problem this has 
already been shown by Paarde kooper & Mellemal (120081) . 

In the following section [2] we describe the physical and nu- 
merical setup of our standard model, and present the numeri- 
cal results in Sect. [3] The validity of the FARGO-algorithm is 
checked in Sect. |4] Alternative setups (nearly local, 2D versus 
3D, adiabatic) are discussed in Sect. [5] The transition of the wake 
into a shock front is discussed in Sect. [6] and in the last section, 
we summarize our results. 



2. The physical setup 

We study planet-disk interaction for planets of very low masses 
that are embedded in a protoplanetary disk. Most of our re- 
sults shown refer to two-dimensional (2D) simulations, using 
the vertically integrated hydrodynamic equations. For validation 
and comparison purposes some additional full three-dimensional 
(3D) models have been performed using a similar physical setup. 
In all cases, we assume that the disk lies in the z — plane and 
use, for the 2D models, a cylindrical coordinate system (r, ip, z), 
while in the 3D case we use spherical polar coordinates (R, <p, 6). 

We consider locally isothermal as well as adiabatic models. 
In the first case, the thermal structure of the disk is kept fixed and 
we chose for the standard model a constant aspect ratio, h = H/r. 
Here r is the distance to the star and H is the local vertical scale 
height of the disk 



H = 



(i) 



where c s is the isothermal sound speed and Ok = (GM./r 3 ) 1 '' 2 is 
the Keplerian angular velocity around the star. During the com- 
putations the orbital elements of the planet remain fixed at their 
initial values, i.e. we assume no gravitational back-reaction of 
the disk on the planet or the star. This allows to formulate the 
problem scale free in dimensionless units. The planet, whose 
mass is specified in term of its mass ratio q = M p /M*, is placed 
on a circular orbit at the distance r — 1 and has angular velocity 
Clp = 1, i.e. one planetary orbit in these units is 27r. The initial 
surface density So is constant and can be chosen arbitrarily as it 
scales out of the equations. 

In the 2D case the basic equations for the flow in the r - <p 
plane are given by equation of continuity 



— + V ■ (Ih) = 0, 
ot 



the momentum equation 
8Ym 



at 



+ V • (lit u) = -VP - SFe 



and the equation of energy 

de 



dt 



+ V • (eu) = -PV u. 



(2) 



(3) 



(4) 



Here, e is the energy density (energy per surface area), and 
P = (y - l)e denotes the vertically integrated pressure. In the 
isothermal case the energy equation is not evolved and the pres- 
sure reduces to P — Zc~, where c s (r) is a given function. The 
external force 

Fext - F* + F p + Finertial- (5) 

contains the gravitational specific forces (accelerations) exerted 
by the star, the planet and the inertial specific forces due to the 
accelerated and rotating coordinate system. 
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Table 1. The physical and numerical parameter for the 2D 
standard model, which consists of a locally isothermal, two- 
dimensional disk with an embedded planet. 



PcLTcimcttir 


Symbol 


Value 


ni occ rnti n 
lllaas laLlU 


n — M I ful 

1 — Mf/m, 


U A 1U 


tlajJCLL IdLlLf 


h = H/r 


U.UJ 


nrvnlinparitv naramptpr 

iiv villi ii v- iiii l y vj ciL 4.1 1 1 1\_ l^i 


M - a [/3 /h 


0.36 


kinematic viscosity 


V 


io- 8 


potential smoothing 




0.1 H 


radial range 


' min ' max 


0.6-1.4 


angular range 


0min 0max 


0-2n 


number of grid-cells 


N r x 


256 x 2004 


spatial resolution 


Ar 


H/16 


damping range at r min 




0.6-0.7 


damping range at r mllx 




1.3-1.4 



For the gravitational force generated by the central star and 
the inertial part of the force we use standard expressions. The 
planetary force is more crucial because it influences for example 
the magnitude of the torque generated by the planet. In our 2D 
standard model we derive it from a smoothed potential and we 
use the following very common form 



T 



2D 



GM„ 



(* 2 + 4f 



(6) 



where s is the distance from the planet to the grid-point un- 
der consideration, and e p is the smoothing length to the oth- 
erwise point mass potential. It is introduced to avoid numeri- 
cal problems at the location of the planet. Then, F p = -VYj; . 
Alternatively, we will use in the 2D simulations a ve rtically av- 
eraged version of F p , following iMuller et al.l ( 1201 2|) . which we 
outline in more detail below. 

Even though we use a non-zero but very small viscosity, we 
do not speci fy those term s explicitly in the above equations, see 
for example iKlevi d!999l) . The viscosity is so small, that it does 
not influence the flow on the relevant scales but is just included 
to enhance stability and smoothness of the flow. In the inviscid 
case, the dynamical evolution of the system is controlled by the 
planet to star mass-ratio q and the pressure scale height H . A 
dimensional analysis by iKorvcanskv & Papaloizoul d 1996b has 
shown that the relevant quantity is the nonlinearity parameter 



M 



,1/3 



2. 1 . The standard model 



The standard model refers to a 2D disk with a locally isothermal 
setup, where the temperature is a given function of radius, that 
does not evolve in time. The parameters for our standard model 
are specified in Table [T] The first two quantities, the mass and 
the aspect ratio, determine the problem physically. As we intend 
to model the linear case in the standard model, we chose a small 
planet with q — 6 x 10~ 6 which refers to a 2ME art h planet for a 
solar mass star. For the disk's thickness we assume h = 0.05. 
For later purposes we list the nonlinearity parameter in the third 
row, which is 0.36 here. All results are a function of q and H 
only, if we assume a vanishing viscosity. In the present situation, 
where we model indeed the inviscid case, we have chosen never- 
theless a small non-zero kinematic viscosity v = 10~ 8 , given in 



units of a p 0^(a p ). This is equivalent to an c-value of 4 x 10 
for a H/r — 0.05 disk, a value that is considered to be much 
smaller than even a purely hydrodynamic viscosity in the disk. 
We have opted for the very small non-zero value for numerical 
purposes. For the planetary gravitational potential ^Fj; we chose 
a value of e = 0.1H in the standard model. We have selected 
this small value for e primarily for c omparison reasons to make 
contact to the recent simulations of iDong et al.l d201 lblal) . who 
suggest a very small smoothing length. Below we will demon- 
strate that for physical reasons a much larger smoothing is re- 
quired, wh ich can be obtaine d by a suitable vertical averaging 
procedure dMiilleret al.ll2012l) . 

A whole annulus is modelled with a radial range from r m [ n 
0.6 to r max = 1 .4. We have chosen this domain size to capture all 
of the torque producing region, which has a radial range of typ- 
ically a few vertical scale heights, see e.g. ID' Angelo & Lubowl 
d2010l) . The computational domain is covered by an equidistant 
grid which has 256 x 2004 grid-cells. This results in a resolution 
of H/16 at the location of the planet for the standard model. To 
reduce or even avoid reflection from the inner and outer bound- 
ary we apply a damping procedure where, within a specified 
radial range, all dynamical variables are damped towards their 
initial values. Specifically , we use the prescription described in 
Ide Val-Borro et al.l d2006l) . and write 



dX 
dt 



X(t)-Xo 



R(r), 



(8) 



where X e {£, u r }, and R(r) is a ramping function increasing 
quadratically from zero to unity (at the actual boundary) within 
the radial damping regions. The relaxation time t is given by a 
fraction of the orbital periods r or b at r m ; n and r max , respectively. 
Here, we use a value of t — 0.03r or b. For more details on the 
procedure see Ide Val-Borro et al 1 (120061) . We note, that it is suf- 
ficient to damp only the radial velocity u r (plus Ug in 3D simu- 
lations) which may be useful in radiative simulations where the 
density stratification may not be know a priori. In test simula- 
tions (not displayed here) that use a larger radial range we have 
found identical results concerning the torque and wake structure 
induced by the planet. 

2.2. Initial setup and boundary conditions 

The initialization of the variables S, T, chosen such that 

without the planet the system would be in an equilibrium state. 
Here, we chose a constant Eo, and a temperature gradient such 
that the aspect ratio h — H/r is constant. That results in T(r) oc 
(7) r~' , which is fixed for the (locally) isothermal models. The radial 
velocity is zero initially, u r = 0, and for u$ we assume a nearly 
Keplerian azimuthal flow, corrected by the pressure gradient 



Jt = 0) = rQ K (l-/ ! / ) 



2x1/2 



(9) 



The planet with mass ratio q is placed at r — 1 and <p - n, i.e. in 
the middle of the computational domain. For some models the 
gravitational potential of the planet is slowly switched on within 
the first 5 orbital periods, while others do not use this ramping 
procedure for the potential. For small mass planets, the results 
that are typically evaluated at 30 r ol -b (i.e. t — 30 ■ 2n), and there 
is no difference between these two options. 

Near the inner and outer radial boundaries the solution is 
damped towards the initial state, using the procedure as de- 
scribed above. In addition we use reflecting boundaries directly 
at rmin and r max . In the azimuthal direction we use periodic 
boundaries. 
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Fig. 1. Density structure of the 2D standard model as generated 
by an embedded planet with q — 6 x 1CT 6 in a disk having the 
aspect ratio H/r — 0.05. Shown is the configuration after 30 T or \,. 
The density is scaled linearly. 




Fig. 2. The total torque, r tot , in units of To (see Eq. [T2l . acting 
on the planet vs. time for the 2D standard model. Shortly after 
insertion, the torque is positive and approximately constant be- 
tween 10 and 40 orbits. At later times it saturates due to mixing 
of the material within the horseshoe region. 



2.3. Numerical methods and codes 

Because one goal of this paper is to verify the accuracy of numer- 
ical methods, we apply several different codes to this physical 
problem. The 2D case is run using the following codes FARGO, 
NIRVANA, RH2D, and PLUTO. All of these are finite volume codes 
utilizing a second order spatial discretization. Additionally, all 
are empowered with the orbital advect ion speed-up, k nown as 
the FARGO-algorithm, as developed bv lMassetl d2000al) . but can 
be used without this algorithm as well. The first 3 codes in th e list 
have been used and described in lde Val-Borro et aL |d200l)The 
last code, PLUTO, is a multi-dimensional Riemann-solver base d 
code for magnetohydrodynamical flows dMignone et al.ll2007l) . 
which has been empo wered recently with the FARGO-algorithm 
dMignone et al.ll2012l) . 

In the standard 2D-setup the simulations are performed in a 
cylindrical coordinate system that co-rotates with the embedded 
planet, and the star is located at the origin. This implies that in- 
ertial forces such as Coriolis and centrifugal force as well as an 
acceleration term to compensate for the motion of the star have 
to be included in the external force term F ext . Additionally, the 
FARGO-algorithm is applied which leads to a speedup of around 
10 for the standard model, and possibly even more for the higher 
resolution models. For testing purposes the FARGO-algorithm 
can be alternatively switched on or off. 

The 3D models are run in spherical polar coordinates. For 
reference we have quoted the inertial terms and their conserva- 
tive treatment in Appendix ICl These are run using the follow- 
ing codes FARG03D, NIRVANA, and PLUTO, where FARGO 3D is a 
newly developed 3D extension of FARGO. For a description see 
AppendixlAl For the timestep we typically use 0.5 of the Courant 
number (CFL). In Sect. |4] we describe the outcome of the com- 
parison in more detail. 



3. Results for the standard case 

To set the stage and illustrate the important physical effects, we 
first present the results of our simulations for the 2D standard 
case using the parameters according to Table Q] For the simu- 
lations in this section we used the RH2D code unless otherwise 



stated. Below, we will then discuss the variations from the stan- 
dard model. 

After the insertion of the planet, the planet's gravitational 
disturbance generates two wakes, in the form of trailing spiral 
arms. The basic structure of the surface density, S, is shown in 
Fig. Q] As seen from the plot, the used damping procedure en- 
sures that reflections by the radial boundaries are minimized. 
There is indication of vortex formation as can be seen by the 
additional structure in the right hand side of Fig. Q] Vortices 
induced by the planet occur for low viscosity disks and have 
been seen already in ear l ier simulations (see e.g. iLi et alj|2005b 
Ide Val-Borro et al.ll2006h In et al.ll2009t) . Here, we will not dis- 
cuss this issue any further. 

The relevant quantity to study the physical consequences of 
the interaction of the embedded planet with the ambient disk is 
the gravitational torque exerted on the planet by the disk. For 
that purpose it is very convenient to calculate the radial torque 
distributi on per unit disk mass, dT(r)/ dm, which we define here, 
following lD'Angelo & Lubowl d2010l) . through the definition of 
the total torque, r tot , acting on the planet 



r* 



In 



r dv 

J dm 



(r) S(r) rdr . 



(10) 



Here, dT{r) is the torque exerted on the planet by a disk annulus 
of width dr located at the radius r and having the mass dm. As 
dT{r)ldm scales with the mass ratio squared and as (H/r)~ 4 , we 
rescale our results accordingly in units of 



dl 



H 



-4 



(ID 



where the index p denotes that the quantities are evaluated at the 
location of the planet, which has the semi-major axis a p . The 
time evolution of the total torque, F tot , is displayed in Fig.|2]for 
the first 500 orbits. The total torque is stated in units of 



H 



(12) 



In this simulation, the planetary potential has been ramped up 
during the first 5 orbits. After insertion of the planet the total 
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Fig. 3. The radial torque density in units of (dT/dm)o (see 
Eq.fTTTi at 30 and 200 r ol -b for the 2D standard model. The torque 
enhancement and spike near r — 1 at t = 30 r or b is due to 
the unsaturated corotation torque. At later times, here shown at 
t = 200 r oi b, only the Lindblad contributions due to the spiral 
arms remain. 
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Fig. 4. Normalized azimuthal density profile of the inner and 
outer wake at radial distances +4/3H away from the planet at 
30 r or b for the isothermal 2D standard model. The coordinates 
x and y refer to local Cartesian coordinates, see Eq. (fT~5T >. The 
'plus' sign in the x-axis label refers to the blue curve at r p -4/3H, 
and the 'minus' sign to the red curve at r p + 4/3//. The upstream 
side of the wake is to the left for both curves. 



torque becomes first positive and remains constant at this level 
for about 30 orbits. In this phase the co-orbital torque, in par- 
ticular the horseshoe drag, is fully unsaturated and gives rise 
to a total positive torque. In our situation of an isothermal disk 
this comes about because of a non-vanishing vortensity gradi- 
ent across the horseshoe region which generates in this case a 
strong positive corotation torque dGoldreich & Tremaindl 19791) . 
The vortensity, which is defined as vorticity divided by surface 
density, is here given by 



(VXH)lz 



(13) 



As can be seen from this definition, for a constant £ disk the ra- 
dial gradient of £ is oc r~ 3 ^ 2 , which leads to the strong vortensity- 
related torque. However, due to the different libration speeds, 
the material within the corotation region mixes and the gra- 
dients of potential vo rticity and entropy are wiped out in the 
absen ce of viscosity (iBalmforth & Korvcanskvl 120011 : iMassej 
1200 ll) . Consequently, the torques drop again, and oscillate on 
timescales on the order of the libration time towards a nega- 
tive equilibrium value, which is given by the Lindblad torques 
generated by the spiral arms. This saturation of the vortensity- 
rel ated to r que re lated torque has been analyzed for example 
by IWardl (120071) through an analyzes of streamlines within 
the horseshoe region for an inviscid disk, a nd later through 
two-dimensional hydrodynamic si mulations bv lMasset & CasolU 
(l2010l) : IPaardekooper et~all(l201 ll) . The strength of the (positive) 
corotation torque depends strongly on the smoothing of the grav- 
itational potential. For the chosen small e = 0.1 this results in a 
positive total torque. For more realistic values of e ~ 0.6 - 0.7, 
r tot will usually be negative, see Section [Sl2l 

To study the spatial origin of the torques we analyze the 
radial torque density for the standard model. In Fig. [3] the 
torque density dT/dm, according to Eq. QT|is displayed vs. ra- 
dius in units of (dT/dm)o. Two snapshots are displayed, one 
at t = 30 r or b where the torque is fully unsaturated and one at 
t = 200 r oi b where the torque is saturated. We note, that for the 
torque calculation we use a tapering function near the planet to 
avoid contributions of material that is bound to the planet, or 



is so close to yield large torque fluctuations due t o numerical 
discre tization effects. We use the form as given in ICrida et ail 
(120081) which reads 



/(*) = 



(14) 



with a tapering length of r t = 0.8 Ru- Here, Rh = (q/3) l ^a p is 
the Hill radius of the planet. Such a tapering is parti cularly use- 
ful fo r massive planets that form a disk around them dCrida et al.l 
120091) . Around lower mass planets, with M < 0.6, circumplane- 
tary disks do not form dMasset et al.l2006l) and a large tapering is 
not required. Indeed, we found that for values of r t in the range 
of 0.4 - 1.0/?h there is not a large difference in the measured 
torques in equilibrium. For example, the variations of the total 
torque in Fig.|2]are less than 5%. 

The torque density in the fully saturated phase, at t — 
200 !T 01 b (Fig. [3] blue line), is positive inside of the planet and 
negative outside of the planet. The positive contribution of this 
Lindblad torque comes from the inner spiral arm, and the nega- 
tive part from the outer one. The distribution at the earlier time, 
t — 30 r oi b, shows an additional contribution and spike just in- 
side of the planet. This part is due to the horseshoe drag which 
is subject to the described saturation process. 

To study the wake properties generated by the planet we 
use here a quasi-Cartesian local coordinate system centered on 
the planet to allow direct comparison to previous linear results. 
Specifically, we define 



x = (r- r p ) and y = {<p- <p p )r p . 



(15) 



In Fig. |4] the relative density perturbations for the inner and 
outer wakes are shown along the azimuth. They are displayed at 
a radial distance of x — +4/3// from the location of the planet. 
For the normalization of the perturbed density we define first the 
thermal mass of the planet 



c 



-H =h 3 M t , 



(16) 
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Fig. 5. The radial torque density of the 2D standard-problem in 
units of (dr/dm)o at 30 r ol -b for different codes at the standard 
resolution. 
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Fig. 6. The radial torque density in units of (dr/dm)o at 30 r ol -b 
for the standard setup for various code. This is enlargement of 
Fig. E] 



where the quantities have to be evaluated at the location of the 
planet. Then, the ratio of the planet mass to the thermal mass is 
given by 

M„ q 



Ma, 



(17) 



Now, we follow lDong et al.l d201 lbl) and scale <5£ = £(</>) - So by 
the planet mass (in units of M tri ) and normalize by x/H. 

Due to the radial temperature variation and the cylindrical 
geometry, the inner and outer wake differ in their appearance. 
However, the general shape and magnitude is very similar to 
the linear result s that have been obtained for the local shearin g 
sheet model (see lGoodman & Rafikovl200lHDong et al.l201 lbl) . 
Differences in the amplitude are presumably due to a different 
normalization. Because our results (in all simulations and with 
all codes) are consi stently by factor of 3/2 larger than those of 
IDong et all d201 lbl) we suspect th at they might have u s ed the 
normalization of M t h as given bv iGoodman & Rafikovl (1200 ll) 
which differs exactly by this factor. At the displayed distance 
from the planet the wake is expected to be in the linear regime, 
which results in a smooth maximum. For this reason we do not 
expect a large dependence on the numerical resolution. In the 
following, we will use only the outer wake to check for possible 
variations due to setup, numerical methods and resolution. 

4. Testing numerics 

To validate our results, and demonstrate that the FARGO- 
algorithm yields accurate results, we varied the numerical setup, 
and used several different codes on the same physical problem. 
In this Appendix we describe our studies in more detail. 

4. 1 . Using different codes on the 2D standard model 

To support our findings on the torque density and wake form 
and demonstrate the accuracy of the used codes, we ran the 2D 
standard model in the isothermal and the adiabatic version using 
all of the above codes. All use the FARGO-setup and are run 
in the same (standard) resolution. The isothermal results for the 
torque density are shown in Figs. [5] and [6] where the latter dis- 
plays an enlargement of the first. Clearly, the results agree ex- 
tremely well for the different codes. This includes the standard 
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Fig. 7. Normalized azimuthal density profile of the outer wake at 
the radius r p +4/3// at 30 T orr) for different codes at the standard 
resolution. 



Lindblad torques as well a the detailed structure of the corota- 
tion torque. The FARG03D code has been used in its 2D version 
for this test. 

Recently, it has been shown by IDong et al] (I201 lbl) : 
iRafikov & Petrovichl (I2012I) that the torque density Y(r) changes 
sign at a certain di stance from the planet in con trast to the stan- 
dard linear results dGoldreich & Tremaindl979l) . Here, we show 
that this effect is reproduced in our simulations as well for all 
codes. In Fig. [6] we show that this reversal occurs at a distance 
r ± « 3.1// a way from the planet in good agreement with the lin- 
ear results o flRafikov & Petrovichl (120121) . Again, all codes agree 
well on this feature, only FARG03D shows small deviations. 

The corresponding wake form at x — r p + AH is displayed in 
Fig. [7] It is identical for all four cases, which shows the consis- 
tency and accuracy of the results and codes. 

4.2. The FARGO treatment 

In Fig.[8]the wake form is analyzed for three different numerical 
usages of the FARGO-algorithm on the 2D standard setup. The 
first, red curve, corresponds to the standard reference case us- 
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Fig. 8. Normalized azimuthal density profile of the outer wake at 
the radius r p + 4/3H at 30 r o ,b for the code RH2D using different 
timesteps and a non-rotating frame. 
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Fig. 9. Normalized azimuthal density profile of the outer wake at 
the radius r p +4/3// at 30 r ol -b using the FARGO-code at different 
resolutions. 



ing a corotating coordinate system and the FARGO-algorithm. 
For the second, blue curve, the simulation has been performed 
in the inertial frame and using FARGO. In the mechanism of the 
algorithm, the quantities in each ring are first shifted according 
to the overall mean angular veloci ty of the ring, and then ad- 
verted using the residual velocity (Masset 2000aj). As a result, 
theoretically it should not matter whether the coordinate system 
is rotating or not. This is exactly what we find in our simula- 
tions, as the blue curve is very similar to the red one. Small 
differences can be produced by the planetary potential which 
is time dependent in the latter case, as the planet is moving, 
and lies at different locations with respect to the numerical grid. 
Also, the number of timesteps used until 30 T OI \, are identical 
for two runs (12, 866 steps). The third, green curve, corresponds 
to a model in the corotating frame without using the FARGO- 
algorithm. Because of the small timestep size in this case, over 
ten times more timesteps had to be used in this case (137,750 
steps). Nevertheless, the wake form is identical. These runs in- 
dicate, that the FARGO-algorithm captures the physics of the 
system correctly. At the same time it comes with a much larger 
timestep and hence much reduced computational cost. This also 
applies to modern Riemann-solvers such as PLUTO, as shown in 
section [5~4l 
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Fig. 10. Normalized azimuthal density profile of the outer wake 
at the radius r p + 4/3// at 30 T OI \, using the FARGO-code with 
different Courant-numbers (CFL). The physical setup differs 
slightly from the standard problem and is described in Sect. 14.41 
In the upper panel the differences of the individual runs with re- 
spect to the standard, CFL=0.5, are displayed. 



4.3. Testing numerical resolution 

In order to estimate the effect of numerical resolution, we ran the 
2D standard model using grid-sizes ranging from 182 x 1418 all 
the way to 4096 x 32064. This is equivalent to grid resolutions of 
H/IQ to Z//256. As shown in Fig. [9] the results are nearly iden- 
tical at all resolutions. The first two, lower resolution cases have 
a slightly lower trough just in front of the wake and a smaller 
amplitude. As will be discussed later in Sect. [6] the results for 
the different resolutions are so similar because at this distance to 
the planet the wake is in the linear regime and has not steepened 
to a shock wave yet. The resolution requirements at the shock 
front will be analyzed in Sect. [6] 

4.4. Testing timestep and stability 

Finally, we would like to comment on possible timestep limita- 
tions due to the gravitational force generated by the planet. In 



our simulations we never found any unsteady evol ution when 
using o rbital advection. In contrast, the results of IDong et al.l 
(l2011bh indicate an unsteady behaviour for larger timesteps. 
They attribute possible instabilities to a violation of an additional 
gravity related timestep criterion and advocate using very small 
timesteps, which would render the FARGO-algorithm not appli- 
cable in very many cases. 

To test specifically this statement we performe d a suite of 
simulat ions on a setup very similar to that used bv IDong et al.l 
(1201 lbl) in their Fig. 12. Due to the difficulty of RH2D and 
FARGO to use a Cartesian local setup we used here a computa- 
tional domain exactly as before with a grid-size of 1024 x 8016 
which gives a resolution of 64 gridcells per scaleheight H. The 
planet mass is 1.33 MEarth* which is equivalent to a mass ratio 
q = 4 x 10~ 6 or M p = 3.2 x 10" 2 M th . For the potential smoothing 
we choose e = 0.0 8// which is yields a planetary potential nearly 
identical to that of lDong et ail d2011bl) . In Fig. [T0]we display the 
results (using FARGO) for different timestep sizes as indicated by 
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Fig. 11. The total torque r tot , in units of To (see Eq. IT2t. acting 
on the planet vs. time. Shown are the 2D standard model (red), 
and a globally isothermal case (blue), with a different density 
profile, such that the potential vorticity gradient vanishes. 

the corresponding Courant-number. The CFL = 0.5 case corre- 
sponds to our standard case. We made the timestep larger (CFL 
= 0.8) as well as smaller, down to CFL = 0.05. All cases yield 
identical results and do not show any sign of instability. In the 
upper panel the differences of the individual runs with respect 
to the standard, CFL=0.5, are displayed. The performed runs 
with the RH2D, FARGO and FARG03D codes yield identical results, 
again with no signs of unsteady behaviour. Here, FARG03D was 
run in the 2D version, both with the setup as indicated above and 
with the local setup of Tab. 3, with resolution /j/64. For all our 
runs, past the first two orbits the wake profile at x — 1.33// has 
achieved convergence to better than the 1 % level, regardless of 
the value of the timestep size. In AppendixlBlwe re-analyze pos- 
sible stability requirements in the presence of gravity and find 
indeed stability for the timestep sizes used with the FARGO- 
algorithm. 

Note, that for these runs we switched off the physical viscos- 
ity completely. We find that the result for the wake displayed in 
Fig.[10]is in fact, due to the special scaling of the axes, identical 
to that of the standard problem as shown in the previous plots. 
Additionally, we have not seen any sign of unsteady behaviour. 
All of this indicates that our small value of the kinematic viscos- 
ity, v = 10~ 8 (in dimensionless units), is essentially negligible. 

5. Using alternative setups 

To illustrate how variations of individual properties of the stan- 
dard model influence the outcome, we performed additional sim- 
ulations which are described in this section. 

5.1. Different radial stratification 

As shown above, in the initial evolution after embedding the 
planet the total torque is positive due to a large positive horse- 
shoe drag. The strength of this effect depends on the radial gra- 
dients of potential vorticity, entr opy (for simulations with en- 
ergy equation), and temperature dBaruteau & Masse! |2012|) . To 
minimize this effect we present an additional, alternative setup 
where the gradients of potential vorticity (vortensity) and tem- 
perature vanish exactly. Hence, we chose a setup with a density 
gradient £ oc r~ 3 ^ 2 and T = const. The time evolution of the to- 



nearly inviscid disks: Numerical treatment 
Table 2. The numerical parameter for the 3D standard model 



Parameter 


Symbol 


Value 


smoothing radius 


' sm 


0.252?h 


radial range 


^min ^max 


0.6 - 1.4 


angular range 


^min '/'max 


0-2/r 


meridional range 


"min ^max 


82°-90° 


number of grid-cells 


N r xN^xNf, 


256 x 2004 x 39 


spatial resolution 


Ar 


HI 16 


damping range at R miu 




0.6-0.7 


damping range at ^ max 




1.3-1.4 



0.1 



0.05 







-0.05 



-0.1 













2D: e=0.1 H 

2D: integrated force 

3D 
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Fig. 12. The radial torque density in units of (dr/dm)o at 30 T m \, 
for 2D and 3D simulations of the standard setup. The red curve 
corresponds to that in Fig. [3] the blue line to a 2D model with 
a vertical integrated gravitational force and the green to the 3D 
model (using NIRVANA). 



tal torque for this model is displayed together with the standard 
case in Fig. QT| Clearly, after the short switch-on period of the 
planet mass, the total torque is negative and constant through- 
out the evolution. This demonstrates that for this density profile, 
Z oc r~ 3 ^ 2 , which resembles (coincidently) the minimum mass 
solar nebula, there is indeed no corotation torque present, and 
the flow settles directly to the Lindblad torque. The final value 
for the total Lindblad torque differs slightly for the two mod- 
els due to the different gradients in density and temperature. We 
note that for this setup, with vanishing vortensity gradient, there 
are also no vortices visible during the initial evolution. 

5.2. Comparing 2D and 3D simulations 

The setup of the described standard case reduces the physical 
planet-disk problem to two dimensions. However, even though 
the disk may be thin, corrections are nevertheless expected due 
its finite thickness. We investigate this by performing full 3D 
simulations using the same physical setup as in the standard 
2D model. The treatment of the inertial forces is outlined in 
Appendix [C] The additional numerical parameters are listed in 
Table|2] The spatial extent and numerical resolution is identical 
to the 2D model. The initialization of the 3D density is chosen 
such that the surface density is constant throughout. In the ver- 
tical direction the density profile is initialized with a Gaussian 
profile as expected for vertically isothermal disks. The temper- 
ature is constant on cylinders. For the gravitational potential of 
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Fig. 13. Normalized azimuthal density profile of the outer wake 
at the radius r p + 4/3H at 30 r ol -b for 2D and 3D simulations of 
the standard setup. The color coding is identical to Fig. [12] 



the planet we chose the so called cubic-form dKlev et al.H2009b 
which is exact outside a smoothing radius r sm and smoothed by 
a cubic polynomial inside of r sm . The advantage of this form lies 
in the fact that in 3D simulations the smoothing is required only 
numerically and the cubic potential allows us to have the exact 
potential outside a specified radius, here r sm . To calculate the 
torque the same tapering function (Eq. [T4b as before has been 
used. 

In Fig.Q~2]we show the normalized torque density dF/dm for 
2D simulations in comparison to a full 3D simulation using the 
same physical setup. Due to the finite vertical extent the torques 
of the 3D model are substantially smalle r than for the corre- 
sponding 2D setup. As lMiiller et al.l (1201 2h have shown recently, 
this discrepancy can be avoided by performing a suitable verti- 
cal averaging procedure of the gravitational force. Specifically, 
the force acting on each disk element in a 2D simulation is cal- 
culated from the projected force that acts in the midplane of the 
disk. Denoting the distance of the disk element to the planet with 
s, the force density (force per area) is given by 



F p O) = - J dz = -GM v s J 



■ dz , 



(18) 



where *T p is the physical 3D potential generated by the planet. 
For the vertical density stratification a Gaussian density profile 
can be assumed for a vertically isothermal disk as a first approx- 
imation. However, the change of the vertical density as induced 
by the planet has to be taken into account. The results using 
this a veraging prescription in an approximate way dMiiller et alj 
1201 2l) is shown additionally in Fig. [12] by the blue curve. The 
overall behavior and magnitude is very similar to the full 3D 
results. For comparison, a 2D model using a fixed e = 0.7H 
(instead of 0.1H of the standard model) for ^Fp yields similar 
amplitude as the full 3D model but a slightly different shape (see 
Fig.[I§. 

In Fig.[13]we show the corresponding wake profile for the 2D 
and 3D setup. For the 3D case, the surface density is obtained by 
integration along the 9 direction at constant spherical radii, R. 
Here, the wake amplitude of the full 3D model is again reduced 
in contrast to the flat 2D case, with e = 0.1. The 2D model using 
the integrated force algorithm yields here again a better agree- 
ment. The 3D results displayed in these plots have been obtained 



Fig. 14. The radial torque density in units of {dF/dm)o at 30 r or b 
for the 3D and 2D simulations of the standard setup. Compared 
are two 3D simulations (using NIRVANA and FARG03D) with a 
2D simulation, using e = 0.7. 



with NIRVANA, but usage of the new code FARG03D yields iden- 
tical results, as is demonstrated in Fig. [14] 

These results demonstrate clearly, that the e-parameter in the 
2D planetary potential pPi 10 ) cannot be chosen arbitrarily small, 
but has to be on the order of the scale height H of the disk. Near 
to the planet a reduction is requi red to account for the reduced 
thickness, see iMiiller et al.ld2012l) . Hence, the value of e = 0. 1 H 
as chosen for the standard setup is too small to yield good agree- 
ment with vertically stratified disks, and serves here only as a 
numerical il lustration to connect to previous linear and nu mer- 
ical results dGoodman & Rafikovll200lt IDong et al.ll201 lbl) . As 
shown bv lMuiler et al.ld2012 ) a value of e — 0.7 H yields similar 
amplitudes to the 3D case, in particular for the Lindblad torque, 
see Fig. [14] However, as can be seen from the figure, for the e- 
potential the relative strengths of the inner an outer torques differ 
from the full 3D and the 2D vertically integrated case. 



5.3. Using a quasi-local setup 

To demonstrate the agreement of o ur simulations with p revi- 
ously published local results, e.g. by IDong et al.l (1201 lblah . we 
have changed the computational setup, which is listed briefly in 
Table [3] Despite the usage of cylindri cal coordin a tes the setup 
is in fact identical to a model used bv lDong et al.l d201 lbl) . The 
very small thickness H of the disk and the small planet mass 
minimize curvature effects and make the problem more local. 
The nonlinearity parameter for this local model is M — 0.32, 
which is similar to the standard case. This quasi-local model has 
been run in a 2D and 3D setup using FARGO 3D. The 3D case 
has been run again in spherical polar coordinates with the same 
spatial resolution as in the 2D setup of Table [3] For the gravi- 
tational smoothing a length of two grid-cells has been chosen, 
which is equivalent here to e = 0.06 H. For the 2D simulations 
we use RH2D and FARGO, while for the 3D simulations we use 
NIRVANA and FARG03D. All these codes are based on the stan- 
dard ZEUS-method and are enhanced with the FARGO-speedup, 
see AppendixlAlfor details. 

In Fig. [15] we compare the torque density of the 2D stan- 
dard model to the quasi-local model. In a local setup any corota- 
tion torques saturate very quickly, possibly due to the very small 
(quasi-periodic) domain in the angular direction. To match this 
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Table 3. The setup for the alternative q uasi-local model. Th e 
parameters have been chosen according to lDong et al.l (1201 lbl) . 

Parameter Symbol Value 



o.i 



-0.05 



-0.1 













2D: standard (Tab. 1) 

2D: local (Tab. 3) 
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Fig. 15. The radial torque density in units of (dr/dm)o. 
Compared is the standard setup with q = 6 x 1CT 6 , h = 0.05 at 
200 r orb to the quasi-local model with q = 3.2x 10~ 8 ,/z = 0.01 
at 30 r or b. The local calculation utilizes the FARGO 3D-code in the 
2D setup. 



condition, the standard model is shown here at 200r or b when the 
corotation torques have nearly saturated. The overall shape and 
magnitude of the two models is qualitatively in very good agree- 
ment, which supports the scaling with (dT/dm)o- For the local 
models a symmetric shape with respect to the location of the 
planet is expected, while for standard model the outer torques 
are larger in magnitude. This can explain some differences. 

In Fig. [16] we compare the wake form of the standard model 
(as shown in Figs. 141 1 3l l to the more local alternative model for 
the 2D setup. The two curves agree very well indeed, despite 
the huge difference in parameters for the planet mass and the 
disk scale height. We attribute the small differences to curvature 
effects. We note that, due the local character of this setup, the 
curves for inner and outer wakes at r p +4/3 look identical for 
the quasi-local model. 

In Fig.[l7]we compare the torque density of the 3D standard 
model to the 3D quasi-local model. As in the 2D case, now the 
overall shape and magnitude of the two models is again quali- 
tatively in good agreement. The local model shows a symmetric 
shape with respect to the location of the planet, as expected. For 
both cases a similar reduction of amplitude in comparison to the 
2D case is seen. 

In Fig. [18] we compare the wake form of the standard model 
to the local alternative model for the full 3D setup. This time 
the two curves for the outer wake agree again very well, despite 
the huge difference in parameters. The profile for the inner wake 
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Fig. 16. Normalized azimuthal density profile of the outer wake 
at the radius r p + 4/3H at 30 r ol -b. Compared is the standard 
setup with q = 6 x 10~ 6 , h = 0.05 to the quasi-local model with 
q = 3.2 xl0~ 8 ,/! = 0.01. 
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Fig. 17. The radial torque density in units of (dT/dm)o for 
full 3D models. Compared is the standard setup with q = 6 x 
10 6 , /i = 0.05 at 30 T olb (using NIRVANA) to the quasi-local 
model with q = 3.2 x 10~ 8 , h = 0.01 at 15 T orb (using FARGO 3D). 

of the standard model deviates from the outer wake as for the 
previous 2D setup. For the local model, inner and outer wake are 
again identical, as expected. 

5.4. Adiabatic simulations 

The assumption of isothermality is only satisfied approximately 
in protoplanetary disks. Because cooling times can be long, it 
maybe more appropriate to take into account the energy equa- 
tion. To study the influence of the equation of state on the out- 
come, we performed purely adiabatic simulations, that solve the 
energy equation (Eq.[4]i together with an ideal equation of state. 
The result of such an approach is presented in Fig. [19] where 
the radial torque density is displayed for the standard isother- 
mal model together with two adiabatic models using y = 1.4 
and 1.01, respectively. The adiabatic results require rescaled 
units because the adiabatic sound speed is by a factor yfy larger 
than the isothermal one. Hence, the pressure scale length is in- 
creased by the same factor, which enters (through H) the units 



mass ratio 


q = 


M p /M» 


3.2 x 1(T 8 


aspect ratio 


h = 
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0.01 


nonlinearity parameter 


M 


= q l ' 3 /h 


0.32 


potential smoothing 
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0.06H 


radial range 
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number of grid-cells 
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Fig. 18. Normalized azimuthal density profile of the outer and 
inner wake at the radii r p + 4/3H for full 3D models. Compared 
is the standard setup with q = 6 x 10~ 6 ,/7, = 0.05 at 30 T olb 
for the outer and inner wake, to the quasi-local model with q = 
3.2 x 10~ 8 , h = 0.01 at 15 T orb , only at the outer wake. 
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Fig. 19. The radial torque density in units of (dF/dm)o (see 
Eq. rXTT > at 30 r ol -b for the 2D standard model for an isothermal 
and adiabatic setup using y = 1.01 and y = 1.40. The units for 
(dT/dm)o and H have been changed for the adiabatic runs, such 
that H -> yH. 
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Fig. 20. The radial torque density in units of (dT/dm)o (see 
Eq. rm > at 500 r or b for the 2D standard model for an isother- 
mal and adiabatic setup using y = 1.01 and y = 1.40. The units 
for (dr/dm)o and H have been rescaled as in Fig. [19] 
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Fig. 21. The radial torque density in units of (dF/dm)o (see 
Eq.nTI) at 30 T mb for the adiabatic standard model using an ideal 
equation of state with y = 1.01. Three different codes have been 
used, RH2D and FARGO are 2nd order upwind schemes and PLUTO 
is a Riemann solver. 



for (dr/dm)o and To. Obviously, there is a huge difference in the 
horseshoe torque between isothermal and adiabatic runs, while 
the Lindblad contributions are similar, once correctly scaled. The 
adiabatic runs yield similar results for the two y values through- 
out. The strong torque enhancement in this adiabatic simulations 
comes from the entropy-related part of the corotation torque 
which is dri ven by a radial gradient of entropy across the horse- 
shoe region dBaruteau & M asset 20081). 

This result is interesting because sometimes an isothermal 
situation is mimicked with an adiabatic simulation using a y- 
value very close to unity. In particular, this may be required by 
Riemann solvers that do not allow to treat isothermal conditions. 
Our results show, that such a n approach has to be treated ver y 
carefully, as shown already bv lPaardekooper & Melle trial (l2008l) . 
They argued that compressional heating near the planet plays an 
important role in determining the torques. Another reason lies in 
the fact that in an adiabatic situation the entropy is conserved 



along streamlines which is not the case for isothermal flows. 
Reducing the value of y even further yields the same results. In 
general, an adiabatic flow with y — > 1 approaches truly isother- 
mal flow only in the the case of a globally constant temperature. 

After a few libration times the horseshoe region is well 
mixed, and the entropy and potential vorticity gradients across 
the horseshoe regions are wiped out. Hence, the horseshoe 
torques disappear and the Lindblad contributions remain. This 
situation is displayed in Figf20]for an evolutionary time of 500 
orbits. Now, the isothermal model agrees well with the adiabatic 
one. 

We have applied several codes on the adiabatic setup as well. 
In Fig.[2T|we display the same results for the adiabatic situation 
using y = 1.01. Again, all codes agree very well, even though 
now the numerical methodology is vastly different as some use 
a second order upwind scheme (RH2D and FARGO) while PLUTO 
uses a Riemann-solver. Only very near to the planet the results 
differ slightly. 
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Fig. 22. The maximum of the density in the spiral wake as a 
function of radius for the 2D isothermal standard model at 30 
T or b. Different numerical resolutions are shown using FARGO. 
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Fig. 23. Normalized azimuthal density profile of the outer wake 
at the radii r v + 5H at 30 T m \, for the 2D isothermal standard 
model. Different numerical resolutions are shown. 



6. Shock formation 

For the damping of the wake it is important where the transition 
to a shock occurs. As a shock indicates a discontinuous change 
in the fluid variables, numerical codes often have difficulty to 
resolve the structure in detail. To analyze this, we plot in Fig. [22] 
the maximum density in the wake as a function of radius for 
various resolutions of the computational grid. At the radius of 
the planet the density obviously has its maximum and it drops on 
both sides. The previous curves for the wake profile have been 
taken near the minimum value of the density maximum. Here, all 
resolutions show an identical maximum of the wake amplitude. 
Hence, as we have demonstrated in Sect. 14.31 the form of the 
wake does not depend very strongly on resolution at a distance 
of \x\ =4/3H. 

Further away from the planet, beyond a distance \x\ > 2H 
the curves begin to differ for the various resolutions. This is 
clearly an indication for non-convergence of the simulations. 
We attribute this to the formation of a shock wave. Indeed, at 
a distance x s ~ 2H from the location of the planet the speed of 
the wake becomes supersonic with respect to the local Kep lerian 
flow. The criterion as given bv lGoodman & Rafikovl(l200lh indi- 



cates for our nonlinearity parameter, M = 0.36, a shock forma- 
tion at a distance of « 2.9H from the planet, consistent with our 
findings. At very high spatial resolutions, i.e. above a grid reso- 
lution of 64 grid-cells per scale height (1024x8016), the curves 
begin to converge in the shock region. Far away from the planet, 
beyond |jc| ~ 6H (r = 1.3) the damping action of the boundary 
condition begins to set in, and the curves coincide again. 

In Fig. [23] the azimuthal density profile is shown at a radial 
locations = r p + 5H- 1.25. At this location the wake is expected 
to have turned into a shock wave. Due to the trailing nature of 
the wake we define the variable y here slightly different as before 
through 

y = (0 - p ) r . 

From the figure it is obvious that the wake has turned into a 
shock at this location. At our standard resolution (256x2004) 
there is no indication for a shock front. The overall form is very 
smooth with the presence of large oscillations behind the wake. 
With increasing numerical resolution the shock becomes better 
and better resolved, but only at the very highest resolution the 
wake turns into a discontinuous jump. The oscillations behind 
the front diminish and move closer to the front with increas- 
ing resolution. Numerical experiments show that these oscilla- 
tions can be damped out by increasing the strength of viscosity. 
However, this smears out the shock front as well. It has been 
suggested that these oscillations are du e to the use d numerical 
scheme and occur for very weak shocks (lReinll2010h . 



7. Summary 

Through a series of 2D and 3D simulations using different com- 
putational methods and codes we have explored in detail the nu- 
merical requirements for studies of the planet-disk problem. In 
our analysis we focus on the torque density acting on the planet 
and the structure of the wake generated by the planet. 

With respect to the applicability of the fast orbital advection 
algorithm, FARGO, we have shown that it leads to consistent 
numerical results that agree extremely well with non-FARGO 
studies. The achievable gain in speed can be significant. For the 
setup used here we found a speed-up of more than a factor of 
10. The method works well in the presence of embedded plan- 
ets, does not show any signs of unsteady behavior, and can be 
applied in two or three spatial dimensions. As it is applicable in 
conjunction with magnetic fields as well, new possibilities with 
respect to numerical s tudies of turbulent accretion disks open up 
dMignone et al.ll2012l) . 

Concerning the treatment of the gravitation al potential of 
embedded planets , we extend previous studies ( lMassetll2003: 
iMtiller et al.ll2012l) to very low mass planets in extremely thin 
disks. We confirm that, for physical reasons, in 2D simulations 
the planetary potential has to be smoothed with about e = 0.6 H 
- 0.7 H. Models where the gravitational force is obtained di- 
rectly through a vertical integration yield always reasonable 
agreement with full 3D simulations. The usage of very small 
smoothing lengths below e = 0.6 H in 2D simulations is not rec- 
ommended, because then the forces in the vicinity of the planet 
are strongly overestimated, which results in an unphysical en- 
hancement of the torque and too strong wakes. 

Through a careful resolution study, we show that the smooth 
wake structure at distances smaller than about 2 H of the planet 
can be resolved well, and consistently, already with very low 
resolution of 8 to 16 cells per scale height. The results are 
clearly converged for 32 grid-cells per H. For larger distances 
from the planet, the spiral wake turns into a shock wave and 
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much higher resolution may be required. We found, that around 
a resolution of about 100 grid-cells per H convergence can be 
achieved. Because this high resolution is only required near the 
spiral shocks, and the flow is relatively smooth outside, numeri- 
cal methods that adaptively refine this crucial region may be the 
method of choice in the future. 

For adiabatic flows we confirm earlier findings 
dPaardekooper & Mellem a 2008) that the unsaturated horseshoe 
drag shows a strong deviation from the isothermal case. Using 
the appropriate scaling the adiabatic corotation torques are 
independent of y and do not converge to the isothermal case, 
even in the limit y — > 1. Hence, the procedure of modelling 
the isothermal case with simulations of y close to unity, has to 
be treated with care. In the final saturated case, where all the 
corotation effects have been wiped out, isothermal and adiabatic 
results agree perfectly, once the correction to the sound speed 
has been applied. 

In Appendix 14.31 we have shown that we do not find an ad- 
ditional timestep criterion due to the planetary potential and we 
also have not noticed any unstable evolution in the case of using 
the orbital advection. The question wh y in the simulation s us- 
ing the ATHENA-code instabilities occur fctong etalj|2011bl) may 
be connected to the trea tment of orbital advection in that code 
dStone & GardinerJl2010h which is apparently different from the 
implementation in the FARGO-code. One should also notice that 
in such simulations the conservative treatment of Coriolis forces 
is ma ndatory to properly conserve angular momentum (IKlevI 
fl998h . 

We have demonstrated that the planet-disk interaction prob- 
lem may be regarded as a very good test to validate an implemen- 
tation of orbital advection, because it admits a nearly analytic 
solution to which a code output can be compared. This is not the 
case for simulations of turbulent disks, where no such known so- 
lutions exist. We hope, that the presented results and comparison 
simulations may serve as a useful reference for other researcher 
in this field. 

Acknowledgements. Tobias Miiller received financial support from the Carl- 
Zeiss-Stiftung. Wilhelm Kley acknowledges the support of the German Research 
Foundation (DFG) through grant KL 650/8-2 within the Collaborative Research 
Group FOR 759: The formation of Planets: The Critical First Growth Phase. 
Some simulations were performed on the bwGRiD cluster in Tubingen, which is 
funded by the Ministry for Education and Research of Germany and the Ministry 
for Science, Research and Arts of the state Baden- Wiirttemberg, and the cluster 
of the Forschergruppe FOR 759 "The Formation of Planets: The Critical First 
Growth Phase" funded by the Deutsche Forschungsgemeinschaft. Pablo Benftez- 
Llambay acknowledges the financial support of CONICET and the computa- 
tional resources provided by IATE. We acknowledge fruitful discussions with 
Ruobing Dong and Roman Rafikov. 



Appendix A: The codes 

For our comparison simulations we utilized the following codes: 
NIRVANA: In its original (F ORTRAN) version a Z EUS -like sec- 
ond order upwind scheme (IZiegler & Yorkd["l997h . with the op- 
tion of fixed nested grids and magneto-hydrodynamics (MHD). 
It can be used in two or three dimensions and can use different 
coordinate systems. Recently, it has been im proved to includ e 
radiative transport and the FARGO-treatment dKlev et al .1120091) . 
RH2D: A two-dimensional radiation hydrodynamics code for dif- 
ferent coordinate systems, origin ally develop ed for treating the 
boundary layer in accretio n disks ( Klevl 19891) . and later adapted 



to the planet disk problem (iKlevll 19991) . 
FARGO: A two-dimensional, special purpose cod e for disk simu - 
lations that first featured the FARGO-algorithm (Masset 2000a). 
The code is publicly available at: http : //f argo . in2p3 . fr/ 



and has been used frequently in planet-disk and related simula- 
tions. 

FARGO 3D: A code based on similar algorithms as the standard 
FARGO-code, but aimed at being more versatile, as it includes 
Cartesian, cylindrical and spherical geometries, in one, two or 
three dimensions, with arbitrary grid limits. Its hydrodynami- 
cal core has been written from scratch, and it includes an MHD 
solver based on the method of characteristics and constrained 
transport. It is parallelized using the Message Passing Interface 
(MPI) and a slab domain decomposition. It is intended in a 
nearby future to run distinctly on clusters of CPUs or GPUs, 
and it will be made publicly available as the successor of the 
FARGO-code. 

PLUTO: A mu lti-dimensional Riem ann-solver based code for 
MHD flows dMignone et al.l 120071) . which can be used in 
the purely hydrodynamic setup as well. Additionally, it 
has been empowere d recently by the FARGO-algorithm 
dMignone et al.l |2012|) . PLUTO is also freely available at: 
http : //plutocode . ph . unito . it 

The first 3 codes in the list have been used and described 
in a n earlier code comparison project on the planet-disk prob- 
lem (Ide Val-Borro et al 1 120061) . There, more massive planets of 
Neptune and Jupiter mass embedded in viscous and inviscid 
disks have been studied for a large number of codes, and the 
focus was on the gap structure of the disk and the total torques 
have been analyzed. 

Appendix B: Timestep limitation in the presence of 
gravity 

Numerically, we expect that possibly gravity might cause prob- 
lems if, due to the gravitational acceleration g, a parcel of ma- 
terial travels more than about half a gridcell of length Ax in one 
timestep At. This requires the additional gravitational criterion 



.Ax 
At G < — 



1/2 



(B.l) 



Using now the smoothed planetary potential of Eq. (|6]l we find 
that the maximum force is given by 



GM n 



(B.2) 



with k = 2/3 3 ^ 2 « 0.4. To obtain the strongest limitation on Af 
we substitute g max in Eq. (IB. lb and obtain 



Ar G < Q K > Wr 1/2 {^pj m . (B.3) 

We compare this limit now to the regular Courant condition 
when using orbital advection which is given by 



At c = — , 



and find 



2 \l/2 



Ar c \kHAx 



(B.4) 



(B.5) 



If there should be no additional timestep limitation generated by 
the gravity then this ratio should be larger than one. Writing now 
for the grid resolution Ax = H/N we obtain finally that 



k H 2 . , 
N>.--M 



(B.6) 
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for stability. With k = 0.4, e = 0.1//, and M = 0.36 we find for 
the necessary resolution N ~ 10. This is indeed fulfilled even 
for our lowest resolution. We point out that this limit formally 
only applies to flows without pressure (dust). If around the planet 
the envelope is hydrostatic, no additional criterion is required. 
Switching on the planetary potential slowly will ensure stability 
throug hout the evolution as will an i nitial atmosphere around the 
planet dDuffell & MacFadvenll2012l) . 

Appendix C: The 3D hydrodynamic equations in a 
rotating frame 

For reference we state here the 3D hydrodynamic equations in a 
rotating coordinate frame. In a coordinate system rotating with 
the (constant) angular velocity Q, omitting pressure terms, any 
external forces (eg. gravitation) and viscosity, the momentum 
equation reads: 

^ + uVu = -2Q x u + i V [(SI x r) 2 ] (C. 1) 

We now use spherical polar coordinates (r, ip, 9), where r is the 
radial coordinate, ip is the azimuthal angle, and 9 is the usual po- 
lar coordinate measured from the z-axis. NOTE, that we use in 
this Appendix the same symbol r for the spherical radial coor- 
dinate. For rotation around the z-axis, SI = Oe ; , the individual 
equations are: 



du 1 

+ uVu r = - [u\ + Mg) + rQ. 2 sin 2 9 + 2u v Sl sin 9 (C.2) 

Urlly UgU ip CO\.8 



dlia, 

— -+UVM„ 

at 



dug U r Ug 
— + uVMfl = + 

at r 



u 2 cot 6 



-2Q. (sin 9u r + cos 9u g ) (C.3) 

+ 2Slu v cos 9 + Qrr sin 9 cos 9 
(C4) 



Introducing the angular velocity oj through 

Uip — r sin 9 lo 
we may write for the three equations dC.2l - [C~4b 

du r 



dt 



+ uVm,. = — + r sin 2 9 (to + Q.) 2 



dUa, 

—— + uVm^ = - (to + 2Q) (sin 9u r + cos 9u g ) 
at 

dug U r Ug , , 

— - + uVug = h r sin 9 cos 9 (to + O)" 

at r 



(C.5) 

(C.6) 

(C7) 
(C8) 



One sees that in the radial and meridional (8) momentum equa- 
tion only the centrifugal part (to + Q.) and in the angular momen- 
tum ((f) equation only the Coriolis term (2Q) occurs. 

C.1. Conservative treatment of Coriolis Terms in Angular 
Momentum equation 



Defining the total specific angular momentum 
h, = r 2 sin 2 0(w + Q) 



(C9) 



and using the continuity equation (in 3D) we may write for the 
angular momentum equation (IC.7b 
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Fig. D.l. Normalized azimuthal density profile of the outer wake 
at the radius r p + 4/3// at 30 r ol -b. Compared is the standard 
setup with q = 6 x 10~ 6 ,/i = 0.05 and the quasi-local model 
with q - 3.2 x 10 8 , / z = 0.0 1 to the linear theoretical results of 
iGoodman & Rafikovl d200lh . 



Expanding h, this may be written: 
d [pr sin 9 (u^ + rCl sin #)] 



dt 



+ V • (ph t VL) = 



(CIO) 



+ V ■ [pr sin 9 (u^ + rQ. sin 6») uj = 

(Cll) 

The validity of this last equation ( IC lib can be easily checked 
by expanding the terms and make use of the continuity equation. 
Then one arrives at the equation dC.7b . In a numerical method 
that evolves u^ the equation ( 1C.1U should be used to solve the 
angular momentum transport conservatively. 

Appendix D: Comparing to linear results 

After submission of the original manuscript, Ruobing Dong 
generously supplied us wit h the data of the linear results of 
IGoodman & Rafikovl (l200lh . In Fig. ID. 1 1 we compare their data 
to our results for the 2D simulations using the standard setup 
of Tab. [T] and the quasi-local setup of Tab. [3] The overall agree- 
ment of our full nonlinear results with the linear case is very 
good. T he small differences between the results are comparable 
to what iDong et aD (1201 lbl) found in their study. Please note, 
that their vertical scaling differs by a factor of 3/2. 
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